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ADAPTIVE ISOGEOMETRIC METHODS WITH 
HIERARCHICAL SPLINES: 

ERROR ESTIMATOR AND CONVERGENCE 

ANN ALISA BUFFA* and CARLOTTA GIANNELLC 


Abstract 

The problem of developing an adaptive isogeometric method (AIGM) 
for solving elliptic second-order partial differential equations with trun¬ 
cated hierarchical B-splines of arbitrary degree and different order of 
continuity is addressed. The adaptivity analysis holds in any space di¬ 
mensions. We consider a simple residual-type error estimator for which 
we provide a posteriori upper and lower bound in terms of local error 
indicators, taking also into account the critical role of oscillations as 
in a standard adaptive finite element setting. The error estimates are 
properly combined with a simple marking strategy to define a sequence 
of admissible locally refined meshes and corresponding approximate so¬ 
lutions. The design of a refine module that preserves the admissibility 
of the hierarchical mesh configuration between two consectutive steps 
of the adaptive loop is presented. The contraction property of the 
quasi-error, given by the sum of the energy error and the scaled error 
estimator, leads to the convergence proof of the AIGM. 


1 Introduction 

The definition of adaptive schemes that provide local mesh refinement is 
an active area of research in the context of isogeometric analysis mm , an 
emerging paradigm for the solution of partial differential equations which 
combines and extends finite element techniques with computer aided design 
(CAD) methods related to spline models. Since the CAD standard for spline 
representation in a multivariate setting relies on tensor-product B-splines, 
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e.g. see [10, 132 ] , an adaptive isogeometric model necessarily requires suitable 
extensions of the B-spline model that give the possibility to relax the rigidity 
of the tensor-product structure by allowing hanging nodes. 

There are a few different frameworks for the definition of splines on 
rectangular tiling with hanging nodes. We mention here T-splines [35, 36j 
that have been used in the context of isogeometric analysis in the pioneer¬ 
ing papers HUE], and their analysis-suitable [33] or dual-compatible mm 
versions. Other possibilities are offered by polynomial splines over (hierar¬ 
chical) T-meshes [T2]H3] or LR-splines [Bid], that have been tested within 
an isogeometric framework in [29] and [22] . respectively. 

Finally, hierarchical splines based on the construction presented in [24] is 
one of the most promising approach. This is also due to the fact that their 
construction and properties are closely related to the ones of hierarchical 
finite elements. Hierarchical B-spline constructions and their use, both as 
an adaptive modeling tool, as well as a framework for isogeometric analysis 
that provides local refinement possibilities, has been recently investigated in 
a number of papers, see e.g. 021 H9l [201 ESJ. 

In the present paper we aim at defining and studying an adaptive isoge¬ 
ometric method (AIGM) based on hierarchical splines. The choice, among 
the adaptive spline models mentioned above, of the hierarchical setting have 
a twofold motivation. On the one hand, it is a natural extension of the 
B-spline model that is able to preserve many key properties directly by 
construction, and the refinement rules are simple and straightforward. In 
addition, although the type of refinement they allow is more restrictive than 
other solutions, the locally structured hierarchical approach allows to defines 
an effective automatically-driven refinement strategy that, in turns, can be 
used to design a fully adaptive method. 

We consider the simple elliptic model problem: 

— div(AV'u) = / in O, u| 9n = 0, (1) 

where H C M rf , d > 1, is a bounded domain with Lipschitz boundary dQ, 
and / is any square integrable function and 

Vx e fl, £ G R d r/i |£| 2 < A(x)£, • £ and |A(x)£| < 772 1 ^| ( 2 ) 

with 0 < r/i < r/ 2 . 

By closely following the framework of adaptive finite elements — see e.g., 
the recent reviews in [301 El] and references therein - for elliptic partial 
differential equations, we aim at designing and analyse the four blocks in 
the following flowchart associated to an AIGM. 
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SOLVE 


ESTIMATE 


MARK 

-> 

REFINE 


At our best knowledge, all previous works on error estimators in isoge¬ 
ometric analysis were mainly devoted to numerical experiments with some 
goal-oriented error estimators based on auxiliary global refinement steps 

mmm- 

Our choices for the different steps of the adaptive loop may be detailed 
as follows. 


SOLVE 


We want to solve problem ([I]) with hierarchical spline spaces. To 
this aim, we define a family of admissible hierarchical meshes, which 
uses the concept of truncated basis [T9], and we consider the Galerkin 
method on these spaces. Admissibility is related to the number of 
levels which are present (with non zero basis functions) on an element, 
and it is a fundamental assumption in our theory. 


ESTIMATE 


We define residual based error estimator for our problem. Thanks to 
the regularity of splines, such an estimator reduces to the L 2 -norm of 
the element-by-element residual suitably weighted with the mesh size. 
We prove that this estimator is reliable, i.e., it is an upper bound for 
the error, and efficient, i.e., it is a lower bound of the error (up to 
oscillations). 


MARK 


We adopt the Dorfler marking strategy m , namely we mark for re¬ 
finement all elements with largest error indicator until a certain fixed 
percentage of the total error indicator is taken into account by the set 
of marked elements. 


REFINE 


A refinement procedure constructs the refined mesh starting from the 
set of marked elements, by following the structure of the recursive 
refine module generally considered in adaptive finite elements, see 
e.g. [271 128], We construct this routine so that the admissibility of 
the refined mesh is preserved between two consecutive iterations of 
the adaptive loop. 


In general, the refinement procedure identifies the mesh with an increased 
level of resolution for the next iteration by rehnining not only the marked 
elements, but also a suitable set of elements in their neighbourhood, anal¬ 
ogously to the concept of refinement patches in an adaptive finite element 
method. This allows to construct a mesh that preserves a certain class of 
admissibility. The refinement mechanism is similar to the strategy adopted 
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to bound the number of hanging nodes per side in the refinement of quadri¬ 
lateral meshes for finite elements [3], and is also related to the properties of 
the domain partitions created by the bisection rule that are needed to prove 
quasi-optimality of adaptive finite element methods PH El M, E9]. 

In the present paper we start the numerical analysis of our AIGM method 
and we provide a convergence result together with the contraction of the 
quasi-error (i.e., the sum of the error and the error indicator), while the 
complexity of the refine routine, together with quasi-interpolation operators 
and optimality of the AIGM, is left to the companion paper [5]. 

The paper is organized as follows. Some preliminary aspects of hier- 
archial tensor-product B-spline constructions are reviewed in Section [2] to¬ 
gether with the definition of truncated hierarchical B-splines (THB-splines) 
and related properties, before introducing the notion of (strictly) admissible 
meshes. The module SOLVE and ESTIMATE of the adaptive isogeometric 
method are discussed in Sections [3] and [4] including an a posteriori error 
analysis in terms of both upper and lower bound for the energy error. Sec¬ 
tion [5] recalls a well-known marking strategy and introduces a refinement 
strategy that preserves the class of admissibility during the iterative loop 
— module MARK and REFINE. Finally, Section [6] concludes the paper by 
summarizing the key results of the present study, and outlines the spirit of 
our companion paper [5], 

2 Hierarchical spline spaces 

We start by considering the hierarchical approach to adaptive mesh refine¬ 
ment, as natural extension of the standard tensor-product B-spline model 
in a general multivariate setting. In particular, we focus on the truncated 
hierachical B-spline basis, since it allows us to identify a certain class of 
admissible mesh configurations. 

2.1 Preliminaries: B-spline hierarchies 

Hierarchical B-spline spaces are constructed by considering a hierarchy of N 

tensor-product d- variate spline spaces V° C V 1 C.C V N ~ l defined on 

a closed hyper-rectangle D in together with a hierarchy of domains P° I) 
P 1 D ... I) P^ -1 , that are closed subsets of D. The depth of the subdomain 
hierarchy is represented by the integer N, and we assume ft N = 0. 

For each level £, with i = 0,1,..., N— 1, the multivariate spline space V f ‘ 
is spanned by the tensor-product B-spline basis B e of degree p = (p\,... ,pd) 
defined on a given tensor-product grid G e . The (non-empty) quadrilateral 
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elements (or cells) Q of G £ are the Cartesian product of d open intervals 
between adjacent grid values. For any coordinate direction i, for i = 1,..., d 
the knot sequences associated to the grids at the different levels contain non¬ 
decreasing real numbers so that each grid value appears in the knot vector 
as many times as specified by a certain multiplicity. At any level £, i.e., for 
the case of standard tensor-product B-splines, the multiplicity of each knot 
may vary between one (single knots) and pi or pi+i for the case of continuous 
and discontinues functions, respectively. In order to guarantee the nested 
nature of the spline spaces V £ C V l+1 , we require that every knot of level 
£ — 1 is also present at level £ at least with the same multiplicity in the 
corresponding coordinate direction. 

From the classical spline theory, it is known that B-splines are locally 
linear independent, they are non-negative, they have local support, and 
form a partition of unity [10,132]. Moreover, there exists a two-scale relation 
between adjacent bases in the hierarchy so that any function s£k ( C V £+l 
can be expressed as 

s = 4 +1 (s)£ ( 3 ) 

/SeS^ 1 

in terms of the coefficients c- 1-1 . 

The domain £1^ is defined as the union of the closure of elements of G £—1 , 
namely 

An element Q of level £ is active if Q C Q £ and any Q* of level £* > £ 
which belongs to any fr +1 ,..., is not a subset of Q. We denote the 

collection of active elements of level £ as 

G e := {Q e G l : Qch 1 A$ Q*e G r , t > £ : Q* C fF A Q* C q} . (4) 

Let Q be the mesh composed by taking the active elements Q at any hier¬ 
archical level, namely 

Q:= {Q£G £ ,£ = 0,...,N-i}. (5) 

For any Q € Q, we define h q := \Q\ l ^ d . A mesh Q* is a refinement of Q if 

each element Q* 6 Q* either also belongs to Q or is obtained by splitting 
Q £ Q in q d elements via “g-adic” refinement, for some integer q > 2. 
The refinement relation between Q and Q* will be indicated as Q* Q. 
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In particular, we will consider the case of standard dyadic refinement with 
<7 = 2. 

A basis for the hierarchical B-spline space can be constructed by a suit¬ 
able selection of active basis functions at different level of details according 
to the following definition, see also [M, T2 \. 

Definition 1. The hierarchical B-spline (HB-spline) basis TL with respect 
to the mesh Q is defined as 

H(Q) := \j3 € B* : supp/3 C A supp/3 % Q* +1 , i = 0, ..., N - l} , 

where supp/3 denotes the intersection of the support of f3 with Q°. 

Remark 2. Note that the hierarchical approach is not confined to dyadic or 
q-adic (uniform) refinement, but it can also handle different kind of mesh 
refinements, including non-uniform configurations. In addition, by assuming 
that the degrees may increase (but not decrease) moving from one level to the 
subsequent in the hierarchy, nested sequence of tensor-product spline spaces 
can be also considered in the context of p- (and k-) refinement. 

2.2 The truncated basis 

We define the truncation of a function s' £ V £ with respect to B c+1 as the 
contributions in © of only basis functions in B £+1 that are passive, i.e., not 
included in the hierarchical B-spline basis I~L{Q). More precisely, 

trunc £+1 s:= ^ c^ 1 (s)/3, (6) 

/3£t3 e+1 , supp /3<£fl e + 1 

where c^ +1 (s) is the coefficient of the function s with respect to the basis 

element (3 at level £ + 1 of the B-spline refinement rule ([3]). By recursively 
applying the truncation to the HB-splines introduced in Definition [lj we can 
construct a different hierarchical basis [El- 

Definition 3. The truncated hierarchical B-spline (THB-spline) basis T 
with respect to the mesh Q is defined as 

T(Q) := {Trunc £+1 0 : fi e B* H H(Q), £ = 0,..., IV — 1 j , 

where Trum/ +1 (3 :=trunc Ar_1 (trunc Ar_2 (... (trunc £+1 (/I))...)), for any f3 € 
B £ nH{Q). 
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The level of a truncated B-spline r G T(Q) is the level of the B-spline 
from which r is derived according to the iterative truncation mechanism 
introduced in Definition [ 3 ] For simplicity, we will denote Tl = 'H(Q), T = 
T(Q ) when there will be no ambiguity in the text. 

2.3 Properties of THB-splines 

The truncated basis T not only spans the same hierarchical space of classical 
HB-splines, namely 

(i) spanT = span H, 

but it also inherits from the hierarchical B-spline basis H. the following 
properties: 

(ii) non-negativity: f >0, Vf eT; 

(iii) linear independence: X^feT c ^ = 0 cy = 0, Vf eT; 

(iv) nested nature of the hierarchical spline spaces for consecutive levels; 

(v) the span of a THB-spline basis defined over a sequence of subdomains 
is contained in the span of a truncated basis defined over a second 
sequence that is the nested enlargment of the original subdomain hi¬ 
erarchy; 

(vi) completeness of the basis: for a certain class of admissible configura¬ 
tions of the hierarchical mesh, span T contains all piecewise polynomial 
functions defined over the underlying grid. 

In addition, the truncation mechanism enriches the THB-spline basis func¬ 
tions so that 

(vii) they preserve the coefficients of the underlying sequence of B-splines; 
(iix) they form a partition of unity; 

(ix) they are strongly stable with respect to the supremum norm, under 
reasonable assumptions on the given knot configuration!]] 

Strong stability of a basis means that the associated stability constants do not depend 
on the number of hierarchical levels. 
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Due the two-scale relation ([3]) between adjacent (non-negative) B-spline 
bases, the non-negativity of truncated basis functions (ii) is preserved by 
construction. Properties (i), (iii)-(v), and (vii)-(ix) are detailed in mm- 
For the analysis of hierarchical spline space in (vi), we refer to [TSJ for the 
bivariate case with single knots, and to [26] for the general multivariate set¬ 
ting with arbitrary knot multiplicities. As a consequence of property (vii), 
quasi-interpolants in hierarchical spline spaces can be easily constructed [33- 
Properties (ii) and (iix) imply the convex hull property, a key attribute for 
geometric modeling applications. 

2.4 Admissible meshes 

The truncation mechanism that characterizes the THB-spline basis can be 
properly exploited to design suitable refinement strategies that define dif¬ 
ferent classes of admissible meshes. A mesh of this kind allows to guarantee 
that the number of basis functions acting on any mesh point is bounded. 
In addition, the support of any basis function acting on a single element of 
an admissible mesh can be compared with the size of the element itself in 
terms of two constants that do not depend on the overall number of hierar¬ 
chical levels. These two properties are the key ingredients for the subsequent 
analysis — see e.g., Theorem [20| related to the a posteriori upper bound, 
and the error indicator reduction provided by Lemma |23| We postpone the 
presentation of the refinement procedure to Section [5j by simply focusing 
here on the desired mesh configuration. 

Definition 4. A mesh Q is admissible of class m if the truncated basis 
functions in T(Q) which take non-zero values over any element Q E Q 
belong to at most m successive levels. 

For this class of admissible meshes, the number of basis functions acting 
on a single mesh element does not depend on the number of levels in the 
hierarchy but only on m, that represents the class of admissibility of the 
mesh. 

Corollary 5. For multivariate tensor-product B-splines of degree p = (pi, 

..., pd), the number of truncated basis functions which are non-zero on each 
element of an admissible mesh is then at most m n?=i(w + i)- 

Another important fact that holds for admissible meshes is the following. 

Corollary 6. If Q is an admissible mesh of class m, given a truncated basis 
function t E T{Q), 

IQI < |supp?| < |Q| VQeQ : Qnsupp?/0, (7) 
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where the hidden constants in the above inequalities depend on m but not on 
t, neither on Q or N. 

In what follows, we will always indicate any inequality which does not 
depend on the depth N of the spline hierarchy with <. 

Since the interplay between the truncated basis funcions that are non¬ 
zero on a certain mesh element and the overall mesh configuration is strictly 
related to the locality of the basis functions, we naturally focus on the support 
extension of an element Q £ . For any fixed level £, the support extension 

collects the elements intersected by the set of B-splines in B e whose support 
overlaps Q. We extend this definition to the hierarchical setting as follow. 

Definition 7. The support extension S(Q,k ) of an element Q E with 
respect to level k, with 0 < k < £, is defined as 

S(Q, k ) := { Q' eG k :3 fie B k , supp fi n Q' ± 0 A supp fi n Q + 0} . 

By a slight abuse of notation, we will also denote by S(Q, k ) the region 
occupied by the closure of elements in S(Q, k). In order to identify a specific 
set of admissible meshes, we also consider the auxiliary subdomains 

Ql '■= U{^ : QcG l AS(Q,£) 

for £ = 0,..., IV — 1. Any u/ represent the biggest subset of Of so that the 
set of B-splines in B e whose support is contained in f fi spans the restriction 
of V e to 

Example 8. A set of admissible meshes of class m = 2 corresponds to 
the restricted hierarchies presented in Appendix A of f and relies on the 
following result. If Of C cD^ 1 for l = 1,..., N — 1, then for any element 
Q £ C) e the THB-splines whose support overlaps Q belong to at most two 
different levels: £ — 1 and £. Figure [7] shows three examples of this class of 
admissible meshes related to the bivariate case of degree {p\ .pi) = (p,p) for 
P = 2,3,4. 

The following proposition generalizes the class of admissible meshes con¬ 
sidered in the previous example to the case of an arbitrary m > 2. 

Proposition 9. Let Q be the mesh of active elements defined according to 
0 and 0 with respect to the domain hierarchy II 0 D fi 1 D ... D Ll N 1 . If 

n*cw*- m+1 , (8) 

for £ = m,m + 1,... ,N — 1, then the mesh Q is admissible of class m. 
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(a) (pi,p 2 ) = (2,2) (b) (pi,p 2 ) = (3,3) (c) (pi,p 2 ) = (4,4) 

Figure 1: Admissible meshes of class 2 considered in Example [8] with respect 
to different degrees. 


Proof. For any t £ T(Q) introduced at level £—m, the function trum/ -?n+1 t 
defined by equation ([b]) is a linear combination of basis functions /3 £ i g f ^ m + 1 
so that = 0. Since 

uj £ -m +1 = \J {<3 : Q £ G l ~ m+1 A S(Q, £ - m + 1) C ^" m+1 } , 

if condition ([8]) holds for £ = m, m+ 1,..., N— 1, then also trumr m+1 = 

0. Consequently, the truncation of a B-spline introduced at level £ — m will 
be non-zero on £l £ ~ m \ This means that any element Q £ Q £ belongs to 
the support of THB-splines of only m different levels: £ — m + !,...,£. □ 

As we will detail later, a relevant set of admissible meshes is the one 
verifying condition ([8]) for £ = m, m + 1,..., N — 1, where different values 
of m > 2 can be considered. 

Definition 10. A mesh Q is strictly admissible of class m if it verifies the 
assumptions of Proposition^ 

The meshes considered in Example [8] are strictly admissible of class 2. 


3 The module SOLVE: the Galerkin method 

In this section we describe our model problem and introduce its discretiza¬ 
tion by means of hierarchical splines. Indeed, we have no aim of generality, 
we consider the most simple elliptic problem. As a first step, we give a 
precise definition of the domain O in which our problem is posed. 
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Given a strictly admissible mesh Qq and the corresponding set of trun¬ 
cated basis function To, we suppose that the computational domain 17 is 
provided as a linear combination of functions in To and control points: 

x E 17, x = F(x) = C<pr(x) x E 17° (9) 

re7o 

where E M d . In all what follows, we suppose that the mapping F : 17° —> 
17 is a bi-Lipschitz homeomorphism: 

||£ Q F|| io0 (no) < C F , ||I7"F _1 || I/ oo( n ) < Cp 1 , |a| < 1 (10) 

where cp and and Cp are independent constants bounded away from infinity. 
We consider then the following problem: 


— div(AVu) = / in 17, u| a ^ = 0, (11) 

where A E C'°°(17) is the diffusion matrix which verifies ©>• 

In order to define the variational formulation of the problem, we consider 
the space of functions in H l { 17) with vanishing trace on <917 

V := -£fo (17) := {v E H\to) : v | 0Q = 0} , 

endowed with the norm \\u\\y = IIVull^^d + IMI^n)- A weak solution of 
@ is a function u E V satisfying 

■u E V : a(u,v) = (/, v), VfEV, (12) 

where a : ¥ x ¥ -> 1 is the bilinear form 


a(u,v) := / AVuVw, Vu, wEV, 

Jn 

and (■,•) stands for the L 2 (17) scalar product. We assume that / E V*. 
The bilinear form a(u, v ) is coercive and continuos with constant ot\ and a. 2 , 
respectively: 


a{u , u) > ai||«||y 
a(u, V) < Q!2||w||v||v||v 


uEV, (13) 

u , v E V. (14) 


Moreover, it induces the energy norm: |||v|||n := a(v, v) 1 / 2 , Vv E V. The 
coercivity and continuity properties of a(u, v ) implies the equivalence be¬ 
tween the energy and the H 1 ( 17) norms on V. In addition, the Lax-Milgram 
theorem ensures the existence and uniqueness of the weak solution (12). 
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We construct now our module SOLVE as the Galerkin discretization of 
(12) by means of hierarchical splines on Q. To this aim, we first need to 


introduce a suitable notation for hierarchical meshes and spaces on Q. 

We consider an admissible mesh Q. such that Q h Qo and we denote by 
T the corresponding basis truncated basis functions. Moreover, we construct 
the corresponding mesh and functions of the physical domain via pullback: 

Q = {Q = F(Q) :Q € Q}. 


For all r E T, we construct: 


r(x) = r(x), x = F(x). (15) 

and we denote by T the collection of all mapped basis functions, and by 
SS(Q ) the space they generate, SS(Q ) = spanT(Q). 

Clearly, Q is a hierarchical mesh on the domain Q and for it, we will 
make use of all the nomenclature introduced in Section[2]by simply removing 
the V First, for all elements Q, we denote by Q its preimage through F, 
i.e., Q = F (Q), and Iiq = \Q\ l ^ d , where \Q\ represents the volume of Q. 
Moreover, we set: 

• Q l = F(Q £ ) and u 1 = F(cD £ ); 

• gt = {Q e q : q e gt} and G l = {Q c fi : Q e G 1 }- 

• for all Q E its support extension with respect to level k is 

S{Q,k) = {Q' e G k : Q'eS{Q,k)}. 


Finally, when Q* is a refinement of Q, we will write Q* V Q, when their 
pre-images Q* and Q verifies Q* V Q. 

We are finally in the position to describe the discrete problem we want 
to solve adaptively. The Galerkin approximation of (12) consists in solving: 

find U E §d(Q) : a(U,V) = (f,V), VVeS d (Q), (16) 


where 

^ D (Q) = {VE55(Q) : V| ao = 0}. 

In the subsequent analysis we assume for simplicity SSd(Q) C C i1 (H). This 
assumption is of course not needed for the development of an adaptive strat¬ 
egy, but it allows us to simplify the analysis, by also showing the specific 
changes with respect to C° finite elements. The general case could be treated 
in a similar way following the classical theory of adaptive finite element 
methods. 
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4 The module ESTIMATE: the residual based er¬ 
ror indicator 

The residual associated to U £ 8 is the functional in V* defined by 

(R,v) := (f,v) - a(U,v), 


that satisfies 


(r,v) = a(u — U,v), Vv£V, 
a(u — U,V) = (r, V) = 0, VEeS. 

By recalling that all discrete functions are continuous with continuous deriva¬ 
tives, we can integrate by parts and obtain 

(r,v)= [ fv — AVC/Vu = f fv - div(AVt/)u, 

Jn Jsi 

where, thanks to our assumption that § C C' 1 (tl), the quantity r = / — 
div(AVCZ) belongs to L 2 (0). In particular, as we expect, this means that the 
residual does not contain any edge contribution as in typical finite element 
indicators El- 

One of the fundamental ingredient in the module ESTIMATE is the 
equivalence between the primal norm of the error and the dual norm of the 
residual: 

I d r > 

||rt — U\\y < —||r||v* < —||u — U\\y. (17) 

aq aq 

As it is standard, in order to use the residual as error indicator, we would 
like to replace the norm || ■ ||y* with the following error indicator 

£q(U,Q) = ^2 £ 2 q (U,Q) with £q(U,Q) = h 2 Q \\r\\ 2 L 2 ( Q y (18) 

QeQ 

When no confusion is possible, we may also abbreviate the above notation 
with e 2 q(U). 

Following [27] (see also El), we will show that the following holds: 

||« - U\W < e Q (U, Q) < ||u - U\W + osc Q (U, Q), (19) 

where 

osc 2 q (U, Q)=J2 osc2 (T, Q ) with osc(17, Q) = h Q \\r - U n r\\ L 2 ( Q) 

Q&Q 
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and n n : L 2 (Q) —» Q n , n = ( 711 , 712 , 713 ), denotes the L 2 projector onto the 
space of polynomials of degree rij in the space direction j. The degrees rij, 
j = 1 ,,d can be fixed large enough so that the oscillation are “smaller” 
than the error |3j. 

Indeed Theorem 12 below, will also provide a local version of the lower 
bound in (19) that reads: 

£q(U,Q) < \\u - U\\y(q) +osc q(U,Q). 

4.1 A posteriori upper bound 

In this section we prove that the residual based error indicator defined in 
(18) is reliable, i.e., it is an upper bound for the Galerkin error. 


Theorem 11. Let u be the exact weak solution of the model problem (12). 
The error of the Galerkin approximation U G §(Q) in (16) is bounded in 
terms of the error indicator £q(U) introduced in (18) as follows: 

\\u — U\\y < C up £q(U), (20) 

where the constant C up is independent on the mesh size and on the level of 
hierarchy. 

Proof. This proof follows exactly the lines of the classical proof of upper 
bound in residual based error estimators. For completeness we repeat here 
the steps that can be found in, e.g., Theorem 6 in m- 


1 


Using (17), we have ||u—I/||y < —||r||y*, and we will prove that IHly* < 

ot 1 

e Q {U). 

Since the basis functions in T form a partition of unity and the residual 
is orthogonal to all basis functions r in T, it holds: 


(r, v) = (r, t v) = ^2 t{v-c t )). 

tgT rer 

By standard Cauchy-Schwarz inequality, we estimate the terms in the 
right hand side as follows: 


{R,t(v-c t )}= f rr(u-c r ) < ||rr 1 / 2 || L 2 (n) ||r 1 / 2 (7;-c T )|| Z/ 2 (f2) . 

Jn 

We denote by ui T = suppT and by h LOr = | supp rf/ d , i.e., its size. We can 
deduce by Poincare inequality that: 

||r 1/2 (u - c r )|| L 2 K) < V-||Vu|| £ a (Wr) d. 
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By taking into account Corollaries [5] and [6j we have 
l - SreT 

2. let h be the piecewise constant function which takes values h(x) = 
\Q\ 1/d , xGQ for all Q € Q. It holds: 


J2 h ^r\\ rTl/2 W ia (wT ) 

r eT 


< 


E 

rer ■ 


>2 2 
h r t 


[ h 2 r 2 = £q(U). 
Jn 


The estimate (20) follows. 


s 


4.2 A posteriori lower bound 


In this section we prove that the residual based error indicator defined in (18) 


is efficient, i.e., it is a lower bound of the Galerkin error up to oscillations. 


Theorem 12. Let u be the exact weak solution of the model problem (12). 


The error of the Galerkin approximation U G S(Q) in (16) bounds the error 


indicator £q(U) introduced in (18) up to oscillations: 


£q{U,Q) < Ci h (\\u-U\\ Y (Q) +oscq(U,Q)), (21) 

where the constant C\ b does not depend on Q. 

Proof. Again, this proof is classical, and we repeat the steps of the proof in 
Theorem 7 of M- First, it is easy to see that 


ll r llv*(Q) ^ l|V(u - U) ||l 2 (q) 
and that the following Poincare estimate is true: 


ll r llv*(Q) < ^qIMIl 2 (<3)- 

Moreover, let <Q) n , n = ( 77 . 1 , 712 ,^ 3 ) be the space of polynomials of degree nj 
in the space direction j, then we know that the inverse inequality holds: 


My* > h Q \\r\\L2( Q) VrGQn, 

where the hidden constant does not depend on Q but it deteriorates with n. 
Finally, if we choose f = II n r, it holds that ||r||v*(Q) < ll r llv*(Q)- 
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Now, all these ingredients can be used together in the following estimate, 
with r = II n r: 


MMU 2 (Q) < h Q \\r\\ L 2 ( Q) + II r - f || L 2 (Q) 

^ lldlv*(Q) + h Q Ik - r\\L*(Q) ( 22 ) 

^ lkllv*(Q) + hQ\\r - f\\ L 2 {Q) 

The proof is completed by setting oscq(U, Q) = hq\\r — r\\ L 2 ^Qy □ 

Remark 13. It should be noted that the lower bound we have proved here 
will not be used in the sequel of the present paper. In fact, contraction of 
the error can be proved without using explicitly the lower bound. We have 
reported this simple proof here in order to collect the main properties of the 
estimator we are using. On the other hand, this lower bound will be needed 
in the companion paper m where optimality will be addressed. 


5 The modules MARK and REFINE 


We now briefly describe the considered marking strategy, before introducing 
a refine module that preserves the class of admissibility of a given strictly 
admissible mesh — see Section 2.4 — and its properties. Finally, we conclude 


this section by discussing the contraction property of our AIGM and its 
convergence. 


5.1 MARK: the marking strategy 

Given an admissible mesh Q, the Galerkin solution U E V(Q), the module 
M = MARK ({e Q (U, Q)} QeQ ,Q), 

selects and marks a set of elements M. C Q according to the so-called Dorfler 
marking [T6] , i.e., by considering a fixed parameter 9 E (0,1] so that 

£ Q (U,M)>9e Q (U,Q). (23) 

This marking strategy simply guarantees that the set A4 of marked elements 
gives a substantial contribution to the total error indicator. 
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5.2 REFINE: the refinement strategy 

The support extension S(Q,k) of an element Q e with respect to level 
k, with 0 < k < £ introduced in Definition [7j analogously, we denote by 
S(Q, k) the support extension of the physical element Q. 

In order to guarantee that a mesh after refinement is admissible, we aim 
at imposing that each active element at level £, Q G belongs to the sup¬ 
port of basis functions of at most levels £ — m + !,...,£. To achieve this, 
given an element Q £ Q e that we want to refine, we select the elements that 
are active, that are of level £ — m + 1 and such that their intersection with 
S(Q,£ — 771 + 2) is not empty. We collect all these elements in a neighbor¬ 
hood of Q, denoted by A f(Q, Q, m) and defined here below in Definition 
Clearly, when Q is refined, all elements in A f(Q, Q, m) have to be refined as 
well and this procedure has to be applied recursively in order to guarantee 
that the final mesh is strictly admissible. 

Definition 14. The neighborood of Q e QnQ e with respect to m is defined 
as 

A f(Q,Q, m) :={q' e £ <-m+1 : 3 Q" € S(Q, £-m + 2),Q" C Q'J , 


14 


when £ — m + 1 > 0, and A7(Q, Q,m) = 0 for £ — m + 1 < 0. 

Figure [2] shows the the neighborood of an element Q with respect to 
777. = 2 when, for simplicity, the identity map is considered. 

An automatic REFINE module which allows to define strictly admissible 
meshes is presented in Figure [3j The core of the refinement strategy relies 
on the internal recursive module REFINE_RECURSIVE. Any element Q on 
which the recursive procedure is called will be subdivided into its children. 
Lemma 15 and Proposition |16| below shows the distinguishing properties of 
this procedure. 


Lemma 15. (Recursive refinement) Let Q be a strictly admissible mesh of 
class 777. The call to Q* = REFINE_RECURSIVE(Q,Q,m) terminates and 
returns a refined mesh Q* with elements that either were already active in 
Q or are obtained by single refinement of an element of Q. 

Proof. For every marked element Q £ Q ( n M. the REFINE_RECURSIVE 
routine is recursively called on any element of level (! = £ — m + 1 that 
belongs to the neighborood of Q with respect to m while £' is greater or 
equal than zero. Since at each recursive call the level £' of interest is 
strictly decreasing, the termination condition will be satisfied after a fi¬ 
nite number of steps. In addition, any element Q touched by a call to 
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1 




(b) (pupi) = (3,3) 







(c) (pi,P 2 ) = (4,4) 

Figure 2: Neighborood A f(Q, Q, 2) (light gray) of an element Q (represented 
by any of the four cells in dark grey) when dyadic refinement is considered 
for some low degree cases. Note that the neighborood is always aligned with 
the grid lines of a previous hierarchical level. 


REFINE_RECURSIVE is subdived in its children only the first time it is 
reached in the return phase after the set of recursive calls. Every element 
of Q is then refined at most once in the refinement process that generates 
Q* from Q. □ 


By exploiting the truncation mechanism in the context of strictly ad¬ 
missible meshes — see Definition 10 — it is possible to show that only the 


supports of truncated basis functions of level i — m + 1,£ — m + 2, 
will contain an element Q E for every refined mesh generated by the 
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Q* = REFINE (Q,M,m) 
for all Q G Q n M 

Q = REFINE_RECURSIVE(Q, Q, m ) 

end 

Q* = Q 


Q = REFINE_RECURSIVE(Q, Q, m) 
for all Q' £ JV(Q, Q , m) 

Q = REFINE JRECURSIVE(Q, Q', m) 

end 

if Q has not been subdivided 
subdivide Q and 

update Q by replacing Q with its children 
end 


Figure 3: The REFINE and REFINE RECURSIVE modules. 


REFINE RECURSIVE module. 


Proposition 16. Let Q be a strictly admissible mesh of class m > 2 and 
let Qm be an active element of level i, for some 0 < l < N — 1. The call 
to Q* = REFINE_RECURSIVE ( Q,QMi m ) returns a strictly admissible 
mesh Q* h Q of class m. 

Proof. Let Ll° D ... D Q N_1 D Ct N , with Q N = 0, be the domain hierar¬ 
chy associated to mesh Q. The refined mesh Q* = REFINETIECURSIVE 
(QiQmP 71 ) contains active elements Q* E G^’* with respect to the domain 
hierarchy fl 0 ’* D ... D D fi N ’* where 

n 0 ’* = n° and I/’* D n l , (24) 


for l = 1,...,IV. Note that the maximum level of refinement in Q* is 


necessarily N according to Lemma 15 

Let Q* E G''’* be an active element of Q*. then Q* C \ for some 

0 < l < N. We have two possibilities: either Q* belongs also to or not. 

• If Q* C then 0 < i < N — 1. Since the initial mesh Q is strictly 

admissible of class m, we have: C u e ~ m+1 , namely Q* C u} e ~ m+1 . 

Now, the refined subdomain hierarchy is a nested enlargement of the 
original one according to (24), and, consequently, w £_m + 1 c 

which implies Q* C 

• If Q* C then Lemma [l5| guarantees that Q* has been obtained 

by applying a single refinement to an element of Q. Hence, there exists 
Q%1 e so that Q%1 — Q* • Condition [s] on Q implies 


Qm c ^ m = 


U{q :Q E G l ~ m f\S(Q,l-m) c H £ - m } 
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and, consequently, 


s(Q%i,e-m) cn E - m 


(25) 


Since is an active element of Q that has been subdivided in 
the refinement process from Q to Q*, the REFINE_RECURSIVE 
module has been called over this element. More precisely, the call 
REFINE_RECURSIVE (Q#, , m) belongs to the chain of recur¬ 
sive calls activated by REFINEJRECURSIVE ( Q,QMi m ) f° r some 
intermediate mesh Q& so that Q* P Q# 8 Q. This mean that the 
recursive routine has been called on any Q' 6 N{Q# , m) with 

= [q' g : 3Q" G S{Q% i ,l-m + 1),Q" C Q'j . 


By combinining (25) with S(Q^,l — m + 1) C — m), we 

obtain S(Q^ /l ,£ — m + 1) C Hence, the coarsest elements in 


S{Q%[,£ — m + 1) are exactly the ones of level £ — m. All these Q’ 


elements of level l — rn have been subdivided into their children of level 
£ — m + 1 in the refinement step from to Q* in order to guarantee 
that 

S{Q* l ,£-m + l)cn < - m+1 -*. 

Then Q* C u e ~ m+1 ’*. 


In both cases, Q* C m+1 ’* implies Q* C m+ 1 >*. Condition Q is then 

satisfied. □ 


The previous results guarantees that the strict class of admissibility of 
the mesh is preserved by the REFINETtECURSIVE module. This result 
extends to the REFINE procedure. 

Corollary 17. Let Q be a strictly admissible mesh of class m > 2 and M 
the set of elements of Q marked for refinement. The call to Q* = REFINE 
( Q,M,m ) terminates and returns a strictly admissible mesh Q* h Q of 
class m. 


Proof. The termination of the REFINE module is directly implied by Lemma 
15 Since every marked element Q activates a call to REFINE_RECURSIVE 
{QiQifn ), in order to prove that the final refined mesh Q* preserves the 
satisfation of ([8]) and, consequently, the class m of admissibility of Q, it is 
sufficient to prove that this property holds after every recursive call. This 
guaranteed by Proposition [l6j □ 


is 


20 




In view of the above corollary, we know that the refine mesh Q* pre¬ 
serves the class of admissibility of the initial mesh Q and, consequently, 
Corollaries [5] and [ 6 ] hold. 


Example 18. An example for the case m 
Figure [^} 


2 and the identity map is shown 




(a) initial mesh (b) marked elements 

Figure 4: The admissible meshes in Fig. [I] are generated by the call to Q* = 
REFINE (Q,M, 2) to the mesh Q depicted in (a) with the marked set M 
of elements shown in (b). 


5.3 Contraction of the quasi-error and convergence 

Following the approach by [ 6 ] (see also im we can prove the contraction 
of the quasi-error , defined as the contribution given by the energy error 
together with the estimator scaled by a positve factor 7 : 

\\\ u ~ ^11 In + 7 £ q(C, Q) 

where, we remind the energy norm is just ||| • |||q = a(-, -) 1 / 2 as defined in 
Section [3j Note that neither the energy error |||it — £/|||n, nor the estimator 
sq(U, Q) considered alone may satisfy a similar contraction property be¬ 
tween two consecutive steps of the adaptive procedure in the general setting 

m- 

In the case of adaptive finite elements, monotonicity of the error is proved 
only under additional assumptions, see e.g., m and [28], and indeed, we 
will not study this property in the present paper. 

We present here only the statement of the contraction theorem. Since its 
proof follows the analogous one for finite elements with only minor changes, 
we postpone it to the Appendix. 
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Theorem 19. Let 6 £ (0,1] be the Dorfler marking parameter introduced 


in (23), and let {Qk, SS(Qk),Uk}k> o be the sequence of strictly admissible 
meshes, hierarchical spline spaces, and discrete solution computed by the 
adaptive procedure for the model problem 0- Then, there exist 7 > 0 and 
0 < a < 1 , independent of k such that for all k > 0 it holds: 

Ilk — £4+111 In + 7 e Q k+1 (^k+h Qk+ 1 ) < « 2 [Ilk — £41 1 In + 7 e Q fc (£k> Qk)\ ■ 

(26) 

An immediate consequence of this theorem, it the convergence of the 
error and of the estimator: 


Corollary 20. Under the same assumption of Theorem 19, both the error 
and the estimator converge geometrically to 0 . I.e., there exists 7 > 0 , 
0 < a < 1 and a constant M such that 

Ilk - £4+i 11 In + 7 e Q fc+ i(£ r fc+i! Qk+ 1 ) < Ma k , 

where M depends on the bounds ([ 2 ]) and (10), but not on k. 

Remark 21. The questions related to convergence for other marking strate¬ 
gies remain open and may require additional assumptions on the refinement 
module which should be further investigated. On the other hand, it seems 
very plausible that any other error estimator verifying the upper bound pro¬ 


vided by Theorem 11_ could be used to replace the simple residual based error 
indicator proposed in this paper. As a side remark, we also note that the 


proof of Theorem [75] does not require the lower bound presented in Theo- 
rem[ 


6 Closure 

A posteriori residual-type estimators for the error associated to the Galerkin 
approximation of a simple model problem have been presented, based on the 
truncated basis for hierarchical splines with respect to some class of admis¬ 
sible meshes and a certain multilevel refinement. In the case of the upper 
bound, two key properties of the basis are exploited together with standard 
inequalities of (adaptive) finite element methods. First, the partition of 
unity property and, second, the bound for the number of basis functions that 
assume non-zero value on any mesh element. In the case of the lower bound, 
classical arguments of finite element estimates can be directly applied. By 
taking into account the a posteriori upper bound previously computed (ES¬ 
TIMATE) and a classical marking strategy (MARK), we introduce a specific 
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refinement procedure (REFINE) to proof the contraction of the quasi-error 
and, consequently, the convergence of the adaptive isogeometric methods. 

Corollary [20] states convergence, but complexity is not analyzed at this 
stage. In other words, we have not proven any connection between the error 
and the number of degrees of freedom that are needed to compute the iterate 
U k . As it is known from the AFEM theory, this can be studied by analyzing 
the complexity of REFINE. In order to do this, we need to understand how 
the refinement module controls the interplay between the number of refined 
elements #Qk — #Qo introduced up to step k (that influences the degrees 
of freedom added during the refinement) and the total number of marked 
elements. Among other things, an estimate of the type: there exists a certain 
constant Ao > 0 such that 

fc-i 

#Qk - #Qo < A 0 £ #Mj 

3=0 

is in need. This kind of complexity estimate has been derived for adaptive 
finite elements in mm for two- and three-dimensional problems, respec¬ 
tively. We will prove an analogous estimate for the adaptive isogeometric 
method here introduced, together with optimal convergence rates, in the 
companion paper [5]. 

Suitable extensions of our adaptive framework may be investigated in 
order to consider less restrictive mesh configurations. For example, the use 
of analysis-suitable T-splines combined with semi-structured hierarchical 
construction has been recently investigated [32117]. The possibility of ex¬ 
tending the adaptivity theory here presented to this case is a challenging 
issue, but the wide modeling capabilities of T-splines encapsulated into the 
hierarchical model would provide a powerful refine module. 
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Appendix 


This appendix is devoted to the proof of Theorem 19 Here we basically 


reproduce the proof of the statement as it was first proved for finite elements. 
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In our presentation, we closely follow Chapter 5 of [31J. We are not adding 
something new, and the value of this appendix is only show that the same 
arguments used for finite elements apply also to our case with very minor 
changes. Indeed, some of the proofs are made easier by the fact that our 
error indicator does not contain jump terms. 


Before proving Theorem 19, we need a few preparatory lemmas. 

This first Lemma is nothing else then the Pytaghoras theorem. See 
Lemma 12 in Ell- 

Lemma 22. Let Q be an admissible mesh and Q* be a refinement of Q, 
i.e., Q* y Q. Let U and U* be the Galerkin solution of problem (16) on 
§£>(Q) and E>d(Q*), respectively. It holds: 


\u-u*\\\l = \\\u-u\\\l-\\\u*-u\\\ 


ci- 


(27) 


Proof. This is an immediate consequence of the Galerkin orthogonality. D 

The next Lemma provides a measure of the reduction of the error indi¬ 
cator £ Q (U) with respect to the mesh Q, see Lemma 13 in [31]. 

Lemma 23. Let Q be a stricty admissible mesh , M. be a set of marked el¬ 
ements and Q* the corresponding refined mesh, i.e., Q* = REFINE(Q, M). 
Then, for allV €S D (Q) it holds MV e V(Q), 

e 2 Q *(U, Q*) < £q(U, Q) - \£q(U,M) (28) 


where 0 < A < 1. 


Proof. For each Q e A4, we denote by Q*(Q) the collection of elements of 
Q* that are created by splitting Q. We know that, by construction, for all 


q* g Q*(Q), it holds hg t < - hfi. Due to (10), there exists then a constant 


c(F), c(F) < 1, independent of Q such that Iiq* < c(F)/ iq. 
If we adopt the notation 


£q*(U,Q)— ^Q*ll r (^)llL 2 (Q*)’ 

Q*eQ*(Q) 


it clearly holds: 


VQeM £q*(U,Q) < c(F) £q(U,Q). 
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Moreover, since the mesh size does not increase for all elements in Q \ A4. 
we have: 

VQ £ Q\M s q *{U,Q) <£q(U,Q). 

Summing up for all Q £ Q, we obtain: 

£q*(U, Q*) < £q(U,Q \ M) + c 2 (F) £q(U,M) 

which implies then p8] ) with A = 1 — c 2 (F). □ 

We turn now to the Lipschitz property of the error indicator sq{U,Q), 
for any Q, with respect to the trial function U. see Lemma 14 in m- 

Lemma 24. Let Q be an admissible mesh and V, W £ § d(Q)- There 
exists a A > 0 such that the following holds for all Q £ Q 

\sq(V,Q) - eq(W,Q)\ < Ar, Q (A,Q)\\V(V -W)\\ L2{Q} , (29) 

where r] Q (A,Q) = /iQ||div(A)|| L oo( Q ) + \\A\\ L oo( Q y 

Proof. By definition, we have: 

r(V) - r(W) = div(AV(F - W)) = div(A) • V(F - W) + A : D 2 {U - W), 

where D 2 - stands for the hessian matrix. Using the inverse inequality 
II D 2 {U - W)\\ L 2 {Q) < h^WVty - W)\\ L 2 iQ) , applying Cauchy-Schwarz and 
triangle inequality, we obtain: 

MV, Q) - e Q (W, Q )I < h Q \\r(V) - r(W) || L 2 (Q) 

< (h Q ||div(A)|| L oc (Q) + ||A|| L oo ( q ) )||V(U - W) || i2(Q) 

which ends the proof. □ 

We can combine the previous results to obtain the last preparatory 
Lemma. See Proposition 3 in }31| . 

Lemma 25. Let Q be a strictly admissible mesh, M be a set of marked ele¬ 
ments and Q* the corresponding refined mesh, i.e., Q* = REFINE(Q, Ai). 
There exists A > 0 so that, \/V £ Sd(Q), V* € § *d(Q*) and any 5 > 0, 

e 2 Q*(V*,Q*) < (1+5) [4(U,Q) - A4(U,M)]+(l+(i- 1 )A 2 7 ? |(A,Q)|||U*-U||| 

(30) 

with 4* =sup Q » eQ * rjQ t (A, Q*). 


25 




Proof. Applying triangle inequality and Lemma [24| we have: 

4* (V*,Q *)<( 1 + 5)4* (V, Q* ) + ( 1 + 5" 1 ) | £q * (V*, Q* ) - £q* (V, Q*)\ 2 
< (1 + 6)eq,(V, Q*) + Q*)A||V(F - V*)\\ 2 L 2( Q *y 

Summing over the elements, we obtain: 


4-(^*> 2*) < (! + 5)4* (V, Q*) + 4*11(1/ - V*)\$. 

The statement follows by applying Lemma [ 


□ 


Finally we are now ready to prove Theorem 19 


Proof of Theorem 19, By summing up the error orthogonality (27) with 


the estimator reduction 
|2 , _2 


scaled by a constant 7 > 0 , we obtain 


u ~ £4+i|||n + 7£ Qfc+1 (£4+i, Qk+i) <|||« — Uk\\\l 


ft 


+ [ 7 (1 + 5 _ 1 )A 0 — l] |||£4_|_i — L/fc|||s 

+ 7 (1 + 5) [£ 2 Qk (U k , Qk ) — Aeg fe (17fc, M. k )\ , 


where we have used (28) with Q = Q k ,Q* = Q k + i,P = U k ,V* = U k +i, 
and we have set Ao = A4 0 (4 Qq) > A rjQ k (A, Qk). The choice 7 = 1/[(1 + 

<f _1 ) Ao] together with the Dorfler marking property (23) leads to 

\\\ u ~ £4+i 11 In + 7 4 fc +i (£4+i) Qk+i) < \ \\u - £4|||n 

+ 7 (1 + 6) [4fc(£4) Qk) — A Q 2 £ 2 Q k (Uk, Qk)\ 

= |||u — £4||ft + 7 (1 + 5)(1 — \0 2 ) e 2 Qk (U k , Q k ). 

By choosing the parameter 5 so that (1 + 5)(1 — A 6 2 ) = 1 — A 0 2 / 2, the above 
inequality reduces to 

( \ 9 2 \ 

1 - 4 fc (£ 4 ) Qk)- 


The second term on the right-hand side may be written as 


A e 2 


-7- 


£ Qi.(Uk, Qk) + 7 ( 1 ~ 


a e 2 


4fc ( U k> Qk), 


so that taking into account the a posteriori upper bound (| 20 |) and the associ- 

1 - 

□ 


ated constant Cup, we obtain the inequality (26) with a = max jl — 7 4 A ^ 2 .1 — 
1 . UP □ 


< 
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